%画出图像的灰度直方图
%画出均衡化后的直方图
%得到均衡化后的图像
clc;clear;
a=imread('pout.tif');
a=a(:,:,1);
subplot(221)
imshow(a);
title('原图');
 b=zeros(1,256);             %构造一个统计灰度出现频次的矩阵
%查找矩阵中各元素出现的次数
for i=0:255              
 b(i+1)= numel(find(a(:,:)==i));
end
c=numel(a);                  %a矩阵中元素总个数
d=zeros(1,256);              %构造一个统计灰度值出现频率的矩阵
%计算a矩阵中元素出现的概率
for j=1:256               
   d(j)= b(j)/c;
end
subplot(222)
bar(d,0.8,'r')               %画出直方图 
title('原灰度直方图')
set(gca,'XTick',(0:10:255))
xlabel('灰阶')
ylabel('频率')    
%灰度直方图二次统计量
m1=0; n2=0;G=0; H=0;
for i=1:256
    m1=(i-1)*d(i)+m1;
    n2=((i-1)^2)*d(i)+n2;
    G=d(i)^2+G;
    H=-d(i)*log2(d(i))+H;
end
fprintf(' 均值%.4f\n 方差%.4f\n 能量%.4f\n 熵%.4f\n',m1,n2,G,H);
%对直方图均衡化
w=zeros(256,8);              %构造矩阵w,1到8列分别是rk,nk,pr(rk),sk计，sk并，sk,nsk,pk(s)
w(:,1)=0:255;                %灰度级
w(:,2)=b;                    %每个灰度出现的次数
w(:,3)=d;                    %每个灰度出现的概率
 w(1,4)=w(1,3);              %先将第一行‘sk计’赋值给’sk并‘
 w(1,5)=round(w(1,4)*255);   %计算第一行变换后的新灰度
%计算变换后的新灰度值 
for ii=2:256
    w(ii,4)=w(ii,3)+w(ii-1,4);  %计算’sk计'
    w(ii,5)=round(w(ii,4)*255);  %计算‘sk并'
end
%计算灰度变换后，每个灰度的出现的频次和概率
for jj=1:255
if w(jj,5)==w(jj+1,5)
    w(jj,6)=0;
    w(jj,7)=0;
    w(jj,8)=0;
    w(jj+1,7)=w(jj,2)+w(jj+1,2)+w(jj,7);
    w(jj+1,2)=w(jj,2)+w(jj+1,2);
else
w(jj,6)=w(jj,5);
w(jj,7)=w(jj,2);
w(jj,8)=w(jj,2)/c;     %计算变换后灰度的频率
end
end
w(256,8)=w(256,2)/c;  %由于循环没有计算最后一行，故需要单独计算最后一行灰度出现的频率
subplot(224)
bar(w(:,8),0.8,'r')             %画出直方图 
title('均衡化后灰度直方图')
set(gca,'XTick',(0:10:255))
xlabel('灰阶')
ylabel('频率')
%画出直方图均衡化后的图像
[m,n]=size(a);
aa=zeros(m,n);
%将变换后的灰度值代替原来的灰度值
for i=1:m
    for j=1:n
        k=a(i,j);
        aa(i,j)=w(k+1,5);
    end
end
aa=uint8 (aa);  %将double型的矩阵aa变成图像的记录格式
subplot(223)
imshow(aa);
title('均衡化后图像');